#include "lencode.hpp"
#include "configs.hpp"
#include <vector>
#include <iostream>

namespace LDPCEncoder{

  // find the index of the first "1" in a vector
  long index_of_first_one(vec_GF2& v){
    long n = v.length();
    for(long i = 1; i <= n; i++){
      if (IsOne(v(i)))
	return i;
    }
    cout << "Bug: the vector contains all zeros, the matrix is not full rank" << endl;
    exit(EXIT_FAILURE);
  }

  // FIXME: this will bring the cache effects on board
  // switch the i-th and j-th columns of matrix m
  void switch_column(mat_GF2& m, long i, long j){
    GF2 temp;
    for(long k = 1; k <= m.NumRows(); k++){
      temp = m(k, i);
      m(k, i) = m(k, j);
      m(k, j) = temp;
    }
  }

  // do column switches if the matrix shape is bad
  void adjust_columns(mat_GF2& H, vector< pair<long, long> >& perm){
    long i;
    long j;
    long m = H.NumRows();
    for(i = 1; i <= m; i++){
      j = index_of_first_one(H(i));
      if (i != j) {
	// start switching columns
	if (DEBUG)
	  cout << "switching columns of H:" << i << " and " << j << endl;
	switch_column(H, i, j);
	perm.push_back(pair<long, long>(i, j));
      }
    }
  }

  // FIXME: works only for a full rank matrix. This function oes inplace
  // conversion to compute the reduced row echelon form of matrix H.
  void to_rref(mat_GF2& H){
    long i, j;
    long m = H.NumRows();
    for(i = 2; i <= m; i++){
      for(j = 1; j <= i - 1; j++){      
	if (IsOne(H(j, i))) {
	  H(j) = H(j) + H(i);
	}
      }
    }
  }

  // get the submatrix from matrix m.
  const mat_GF2 submatrix(mat_GF2& m, long i1, long i2, long j1, long j2){
    mat_GF2 result = mat_GF2();
    result.SetDims(i2-i1+1, j2-j1+1);
    long a, b, p, q;
    for(a = i1, p = 1; a <= i2; a++, p++){
      for(b = j1, q = 1; b <= j2; b++, q++){
	result(p, q) = m(a, b);
      }
    }
    return result;
  }

  // compute the generator matrix from the parity check matrix
  void compute_G(mat_GF2& G, mat_GF2& H){
    long n = H.NumCols();
    long m = H.NumRows();
    long k = n - m;

    G.SetDims(k, n);
    mat_GF2 Hcopy = mat_GF2(H);
    vector< pair<long, long> > perms;
    gauss(Hcopy);
    adjust_columns(Hcopy, perms);  
    to_rref(Hcopy);
    mat_GF2 transposedMat = transpose(submatrix(Hcopy, 1, m, n - k + 1, n));
    long i, j;
    for(i = 1; i <= k; i++){
      for(j = 1; j <= m; j++){
	G(i, j) = transposedMat(i, j);
      }
    }
    // assign identity diagonal matrix to G
    for(i = 1; i <= k; i++){
      G(i, m + i) = to_GF2(1);
    }
    // perform inverse column permutation
    for(i = perms.size() - 1; i >= 0; i--){
      switch_column(G, perms[i].first, perms[i].second);
    }
  }

  // encode message with generator matrix G.
  void encode(vec_GF2& codeword, const vec_GF2& msg, const mat_GF2& G){
    codeword = msg * G;
  }


  void decode(vec_GF2& msg, const vec_GF2& codeword, const mat_GF2& G){
    // to be implemented, in this work this function is not necessary
  }


  // generate a random codeword
  void generate_random_codeword(vec_GF2& codeword, const mat_GF2& G){
    codeword.SetLength(G.NumCols());
    srand(time(NULL) * rand());
    sleep(0.2);
    SetSeed( to_ZZ(rand()) );
    encode(codeword, random_vec_GF2(G.NumRows()), G);
  }

  // verify the correctness of the generator matrix G
  void verify_G(const mat_GF2& G, const mat_GF2& H){
    vec_GF2 codeword;
    generate_random_codeword(codeword, G);
    vec_GF2 x = H * codeword;
    if (IsZero(x)){
#ifdef VERBOSE
      cout << "passed verification" << endl;
#endif
    }
    else{
      cout << "bug in G matrix gengeration" << endl;
      cout << x << endl;
      exit(EXIT_FAILURE);
    }
  }

} //LDPCEncoder
